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Abstract. The Nonlinear Schrodinger Equation (NLSE) with a random 
potential is motivated by experiments in optics and in atom optics and is a 
paradigm for the competition between the randomness and nonlinearity. The 
analysis of the NLSE with a random (Anderson like) potential has been done 
at various levels of control: numerical, analytical and rigorous. Yet, this model 
equation presents us with a highly inconclusive and often contradictory picture. 
We will describe the main recent results obtained in this field and propose a list of 
specific problems to focus on, that we hope will enable to resolve these outstanding 
questions. 
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1. Introduction 

The Nonlinear Schrodinger Equation (NLSE) with a random potential is a 
fundamental problem. In spite of extensive mathematically rigorous, analytical and 
numerical explorations, the elementary properties of its dynamics are not known. The 
problem is relevant for experiments and its resolution will shed light on many problems 
in chaos and nonlinear physics. It may also stimulate novel experiments. On a one- 
dimensional lattice, the NLSE with a random potential (that will be the subject of 
the present review) is given by, 

id t tp{x,t) =H iP(x,t) + f3\ib(x,t)\ 2 ib(x 1 t), (1) 

where 

H ^p {x, t) = -J [V» (ar + 1, t) + tp (x - 1, t)} + e x ^ (x, t) , (2) 

while, i£Z; and {s x } is a collection of i.i.d. random variables uniformly distributed 
in the interval [— ^, ^ ]. The Hamiltonian Hq is the Anderson model in one-dimension 
[TJ [21 [31 131 [5] . It is important to note that for the dynamics generated by ([1) the £ 2 
norm, J\f — J2 X IV* i x i ^)| 2 ' an d the energy given by the Hamiltonian ([8]) are conserved 

In the present review we will focus on the question about the dynamics; will an 
initial wave function ij)(x, t = 0), which is localized in space, spread indefinitely for 
large times, and in particular in the asymptotic limit, t — > oo. Surprisingly, the answer 
to this elementary question is not known in spite of extensive research in the last two 
decades [3 [HI El [TOJ [HI [HI H31 H3] ■ We believe that the NLSE © is a representative of 
many nonlinear problems, such as the famous Fermi-Pasta-Ulam (FPU) problem |15| . 
Therefore, its understanding may shed light on dynamics generated by other nonlinear 
equations, e.g, nonlinear Klein-Gordon and FPU equations. Many properties of (fTJ) 
will be shared by the continuous counterpart, where a; is a continuous variable. Since 
most of the results that were derived so far are for the discrete problem, in the present 
review we will confine ourselves to this case. The dynamics is completely understood 
in the two limiting cases. In the absence of the random potential (e x — 0, for all 
x) an initially localized wavepacket will spread indefinitely, for all values of (3, unless 
solitons are formed. In the discrete case, unlike the continuous case, the formation of 
the solitons cannot be established rigorously [16]. The continuous version of this model 
is in fact an integrable problem [6J. For attractive nonlinearity, ft < 0, solitons are 
found while for repulsive nonlinearity, /3 > 0, complete spreading takes place. In the 
presence of randomness (w > 0), but for (3 — Eq. ([T]) reduces to the Anderson model 
© where it is rigorously established that all the eigenstates are exponentially localized 
in one-dimension with probability one [1] [2j [4j [5] . At long scales the eigenfunctions 
behave as 

u n (x) ~ e -l— (3) 

where x n is the localization center and £ is the localization length. Consequently, 
diffusion is suppressed and in particular a wavepacket that is initially localized will 
not spread to infinity. This is the phenomenon of Anderson localization. In two- 
dimensions it is known heuristically from the scaling theory of localization, that all the 
states are localized, while in higher dimensions there is a mobility edge that separates 
localized and extended states [3J[3]. 

The behavior of the dynamics generated by (flj is very different in the two 
extreme limits (w = 0, ^ 0) and (w ^ 0, ft = 0). Therefore, it is a paradigm for the 
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exploration of the competition between randomness and nonlinearity. Let us comment 
on the choice of the nonlinearity. 

The nonlinear term of the form /3 \ip\ 2 ip used in ([1} is just one possibility, which 
is used in this review for the sake of clarity. In several theoretical studies it is replaced 

by nana ma hbi 

H a =(3\ij\ 2 ^, (4) 

where a > is arbitrary. In some mathematical studies it is also replaced by 
Hp(x) — P {x) \ip\ 2 ipi where j3 (x) is a decaying function of x [7J. Other types types of 
nonlinear terms appear in experimental realizations. 

The NLSE was derived for a variety of physical systems under some 
approximations. It was derived in classical optics, where ip is the electric field, by 
expanding the index of refraction in powers of the electric field, keeping only the 
leading nonlinear term [15]. Let the index of refraction depend on the intensity of the 
electric field \ip\ , then for weak fields it takes the form, 

n(|^| 2 ) =n + m|V| 2 + o(M 4 ) ■ (5) 
The nonlinear term in (JTJ corresponds to a weak field so that the quartic 
correction is negligible. In several important cases n (\ip\ 2 ^ saturates, namely, 

lim^p^^ n (\ip\ 2 ^J — const. For example in the induction technique [201 121] the 
index of refraction takes the form, 

n 



2- (6) 

In optics, Eq. (JTJ) corresponds to the paraxial approximation where the 
propagation direction plays the role of time. In this approximation the variation 
in the index of refraction in space is weak, and therefore there is only a small change 
in the propagation direction, and back-scattering is negligible. 

For Bose-Einstein Condensates (BEC), the NLSE is a mean field approximation, 
where the term proportional to the density /3|i/'| 2 approximates the interaction between 
the atoms. In this field the NLSE is known as the Gross-Pitaevskii Equation (GPE) 
(25] H31 122 I2FJ I2S1 [27] • It was rigorously established, for a large variety of interactions 
and of physical conditions, that the NLSE (or the GPE) is exact in the thermodynamic 
limit [28, 29J. Experiments on spreading of wavepackets of cold atoms in a random 
optical potential were recently performed |30[ I5T1 [521 155] . In those experiments as 
in experiments in optics, the random potential exhibits correlations and therefore 
deviates from the model presented in ([T]). 

Another possible form of the nonlinear term is 

H H =ip(x)Jv{x-x')\iP{x')\ 2 , (7) 

which results in the Hartree approximation extensively used in solid-state physics. The 
Gross-Pitaevskii equation (or NLSE) is obtained from ([7]) in the limit where \ip (x')\ 2 
is slowly varying. 

The theory of Anderson localization was very recently extended to the many-body 
particle systems [341 [551 [55] , It was found that indeed for sufficiently low energies 
localization takes place for fermions |35[ [55] as well as for bosons [37]. It should 
be emphasized that in these works localization is analyzed for the case where the 
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density is non- vanishing in the thermodynamic limit. The problem of spreading with 
a vanishing average density, corresponding to the problem that is the subject of the 
present review, is different in principle in the thermodynamic limit. For the N body 
problem, a wavepacket that is initially localized will remain localized, as established 
rigorously in [38, 35] (see Subsection 13. ip . 

The model ([1]) was motivated by experimental realizations we have discussed 
above, but this review will treat the problem of spreading, and in particular the 
asymptotic one, as a fundamental theoretical problem. 

For linear problems, all aspects of dynamics are determined by the spectral 
properties, namely the eigenvalues and the eigenfunctions. This is not correct for 
nonlinear problems. For example, for small /3 in (QJ there are stationary and quasi- 
periodic states which are exponentially localized |4l)[ |4"TI 1^2"] . This however does not 
imply that an initially localized wavepacket will not spread, contrary to the case of 
a linear system with a bounded localization length. Transmission through a chain 
fl} was extensively studied [33J S31 ESI US] , but since it is not directly related to the 
spreading problem (unlike the situation for linear systems), it will be left out of the 
scope of the present review. 

The natural question we will survey is whether a wavepacket, that is initially 
localized in space, will indefinitely spread for dynamics controlled by ([1]). A simple 
heuristic argument indicates that spreading will be suppressed by randomness. If 
unlimited spreading takes place, the amplitude of the wave function will decay since 
the £ 2 norm, A/", is conserved. Consequently, the nonlinear term will eventually become 
negligible, and Anderson localization will take place as a result of the randomness, as 
was conjectured by Frohlich et al |47| . However, in numerical calculations performed 
by Shepelyansky [3B] for the kicked rotor with a cubic nonlinear term, Anderson 
localization (that takes place in the absence of the nonlinear term) was destroyed 
and sub-diffusion takes place. Similar spreading was found numerically also by 
Shepelyansky and Pikovsky [8] and by Flach and coworkers \12\ 110). Therefore, the 
naive argument for localization of ([1]) has to be reconsidered and a proper theory 
should be developed. A natural question is what can we conclude from the numerical 
simulations ? The main problem is that dynamics of ([!} are chaotic. The dynamics 
are generated by the Hamitonian, 

#nls (V (x, *)) = E t J ^ ( x + !) ^* (*) +i>*{x + l)iP (x)) (8) 

X 

+ e x \^{x)\ 2 + ^{x)\ A 

Where the NLSE ([T]) is the corresponding Hamilton's equation with the conjugate 
variables ip (x) and ip* (x). Due to the nonlinearity, the motion in the ip(x),ip* (x) 
phase-space will be typically chaotic. Therefore, the numerical solutions of ([1} are 
not the actual solutions. In order to draw conclusions it is assumed that they are 
statistically similar to the correct solutions. Since it is a system of an infinite number 
of degrees of freedom there is no real theoretical support for this assumption. If 
we use the fact that only a finite number of the ip (x) variables are involved, there 
is a competition between two effects. Chaos is enhanced by increase in effective 
number of degrees of freedom, and suppressed by the decreasing amplitude of the 
spreading wavepacket. This competition is outlined in Subsection 12.21 There may be 
also technical problems with the numerical algorithm, this will be discussed in Section 

in 
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Another model similar to |T]), that was extensively studied, is the quartic Klein- 
Gordon equation [T^l [TU] , 

#KG = V \p 2 x + \e x q 2 x + -q% + - (q x+1 - q x f . (9) 
z — ' 11 4 w 

X 

This model differs from the NLSE with a random potential, where the modes of 
the linear problem are effectively localized on a few lattice sites. The perturbation 
theory is well controlled and Nekhoroshev type estimates were established [49] ■ 

A central problem in the field is the Fermi- Pasta-Ulam (FPU) problem [T5]. It 
is interesting to point out that for this problem the spreading starts after a long 
time, although it is very different from the NLSE and the nonlinear Klein-Gordon 
problems. For the NLSE and the quartic Klein-Gordon equations, only neighboring 
modes in space, of the corresponding linear problem are coupled, while for the FPU 
problem, all the modes are coupled. 

2. Theoretical analysis 

In this section various non rigorous theories for the spreading mechanism will be 
presented. The purpose is to analyze the various regimes starting from the short time 
regime and up to the asymptotic time regime. Various authors use an expansion of 
the wavefunction in terms of the eigenstates, u m {x), and eigenvalues, E m , of Hq as, 

i } { X ,t)=Y J Cm{t)^ iEmt Um{x). (10) 
m 

For the nonlinear equation the dependence of the expansion coefficients, c n (t) , is 
found by inserting this expansion into (J]), resulting in 

id tCn = p v; mim2m3 c ; ilCm2 c m3e < fJ " +fJ -^ fi -- B '"3)* (ii) 

mi ,T7l2 ,7713 

where V^ nirn2rrl3 is an overlap sum 

X 

This sum is negligibly small if the various eigenfunctions are not localized in the same 
region of the order of the localization length, £. 

The eigenvalues E n , the eigenfunctions^n (x) , the expansion coefficients, c n (t) 
and the overlap sums depend on the random potentials, {e x } and therefore they 
are random varibales. Consequently, E n , u n {x) and V™ 1 ™ 2 ™ 13 take different values 
for the various realizations of the potentials, {e^}. In Subsections 12.11 and 12.21 
phenomenological theories are discussed, and in Subsection 12.31 a relation to phase- 
space structures is briefly presented, while in Subsection 12 . 41 a systematic perturbation 
theory is developed. 

2.1. Effective noise theories 

In this subsection we present the phenomenological theory |121 [TU] for the spreading 
that is found numerically and will be described in detail in Section @] It is clear 
that not all the content of the initial wavepacket spreads [50]. It was rigorously 
shown, that for sufficiently large /3, the initial wavepacket cannot spread so that its 
amplitude everywhere vanishes at infinite time |50) . It does not contradict spreading 
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Figure 1. (color online) Probability distribution |i/>| over lattice sites x for 
J = 1, w = 4 and for P = 1, t = 10 s (top blue/solid curve); t = 10 5 (middle 
red/grey curve); = 0, t = 10 5 (bottom black curve). (Fig. 2 of l8l). 



of a fraction of the wavepacket. For a wavepacket initially localized on one lattice 
site or started at one linear eigenstate of H , sub-diffusion was found in numerical 
experiments [48, 5TJ [501 EH H2] • The purpose of the theory presented in what follows, 
is to explain the spreading that takes place after some time. It was found numerically 
that after some initial time the shape of the wavepacket is similar to the one found in 
Figffl 

It consists of a relatively flat region at the center and exponentially decaying 
tails. The theory of [TH [23 assumes spreading from the relatively flat region of the 
wavepacket to the region where the amplitude of the wavepacket is small. Let mi, ?7i2 
and TO3 designate eigenstates of Ho with the centers of localization found within the 
flat region, and let n designate a state with a center of localization found in the tail of 
the wavepacket, but in the vicinity of the flat region. Therefore, spreading will take 
place to the region where the n-th state is localized, 

\c mi I 2 ~ |cm 2 1 2 ~ |c m3 1 2 w p (13) 

while 

|c„| 2 <p. (14) 

It is further assumed that the RHS of ([TTj) is a random function denoted by F n (t) , of 
the form 1121. 



where 



F n = dPPp^Ut) 



P = C 2 (3p. 



(15) 



(16) 



Equation (116)) is the probability of a "resonance" between four modes, C\ and C2 
are constants, while /„ (t) is a random function with a rapidly decaying correlation 
function, C (t) . Introduction of P and the assumption that it satisfies (fT5|) , in 
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particular its linearity in /3, are the strongest assumptions of the theory, which still 
require justification. In [TJ] it is claimed that numerical calculations support this 
assumption. Under these assumptions (ITT|) reduces to 

id t c n {t) =F n {t). (17) 

Assuming that F n (t) can be considered random, with rapidly decaying correlations in 
time, integration results in 

c„(i) = -idPPp^ 2 f dt'f n (t'). (18) 



Averaging over realizations one finds 

(|c„| 2 ) =C 3 p* P H, (19) 

where 

C 3 =CxC 2 / C(t') dt'. (20) 
Jo 

The equilibration time, T, is the time when (icnj 2 



T «^4- (21) 



1 

During the time, T, the flat region spreads so that it covers the site where the state n 
was localized. The resulting diffusion coefficient satisfies, 

D oc ^ oc /?V- (22) 

At time scales t 3> T, diffusion takes place and 

^- oc M 2 oc Dt oc ^pH (23) 
p 2 

Where Mi — ^ x \ip (x, t)\ 2 is the first moment and 

M 2 = ^(a;-Mi) 2 |V(x,t)| 2 (24) 

X 

is the second moment. Therefore, 

1 ~ ^4 A 1/3 



->2 



3C 



(/3 4 i) 7 , (25) 



P 

and the second moment satisfies, 

M 2 oc /3 4 / 3 ^/ 3 . (26) 

In agreement with the numerical results presented in Fig. [5J 
The equilibration time satisfies, 

T oc r i/3 t 2/3 . (27) 

Therefore the theory is consistent, as T <C t for large t. 

The crucial assumption of the theory presented in this section is that F n (t) 
behaves as noise with a rapidly decaying correlation function so that the integral 
(|20p converges. This assumption was explicitly tested [52J. The reasons for F n (t) to 
behave as a random variable is that the sum (|11[) consists of many terms with random 
phases, and the dynamics of the c n (t) are chaotic, since these are generated by the 
nonlinear Hamitonian, H^^g. 

A crossover to the regime (|26p from the regime where a different power law is 
found is presented in [53"] . 



The Nonlinear 




Figure 2. (color online) M2 versus time in log-log plots. For J = 1, w = 4 and 
/3 = 0,0.1,1,4.5 ((o)range, (b)lue, (g)reen, (r)ed). The disorder realization is kept 
unchanged. The dashed straight line guides the eye for exponent 1/3 (Fig. 2 of 

El). 



2.2. Scaling 

The theory presented in the previous subsection assumes that F n (t) is effectively 
random inside the relatively flat region of the wavepacket. Chaotic behavior of the 
ip (x) would generate such an F n . As time evolves, the width of this region increases, 
and the density, p, decreases. The first is expected to enhance chaos, and the second 
to suppress it. The calculation of the present subsection [S3] was performed in order to 
check which effect wins. Since one is limited in the length of the system one can study 
numerically, the calculation is performed for a system that is finite but its behavior 
was found not to change as the size of the system is increased. Therefore, it can 
be extrapolated to arbitrarily large system sizes. In particular, the probability for a 
system to be regular as a function of its length, L, density, p and the distribution 
of the random potential, in, was studied. A system is considered regular, if all its 
Lyapunov exponents vanish, and it is considered chaotic, if at least one Lyapunov 
exponent is positive. In the present subsection a theory for the flat region of Fig. |T| is 
presented. For this purpose a dynamics generated by ([1]) for a finite system of size L 
is examined. The linear system is invariant under the rescaling, 

J -> J' = cJ, w ->w' = cw (28) 

if time and energy are correspondingly rescaled. In particular, the localization length is 
unaffected. For the nonlinear systems also the rescaling of the norm or the nonlinearity 
coefficient, 

/3 -> = c/3, (29) 
is required. In the present work the choice 

C -7^TW) (30) 
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was made resulting in, 

J' = — !— w' = — — P' = — ^ -. (31) 

1 + W J (1 + W) H J (1+W) K ' 

In particular we may choose J = 1, and W = w/2 and study the model (JTJ with J 
replaced by J' = 1/ (1 + W) the random potential satisfies, 

W W 

<e x < -fl— (32) 



1 + W ~ 1 + W 
and j3 replaced by f3' = f3/ (1 + W). For convenience this factor can be absorbed in 
the definition of the norm of the wavefunction, 

M ' = TVw' 

The density is, 

It was found that for the fixed W and p, the probability to observe regularity decreases 
with the length L. This behavior can be understood as follows. Suppose we fix W, p 
and consider a lattice of large length L. One divides this lattice into (still large) 
subsystems of lengths Lq. How the probability to observe regularity on the large 
lattice P(p, W, L) is related to the corresponding probabilities for smaller lattices 
P(p, W, Lq)? It is reasonable to assume that to observe regularity in the whole lattice 
all the subsystems have to be regular, because any chaotic subsystem will destroy 
regularity. This immediately leads the relation 

P( P ,W,L) = [P(p,W,L )] L / L ° . (35) 

Equation (|35p implicitly assumes that chaos appears not due to an interaction between 
the subsystems, but in each subsystem (of length Lq) separately. This appears 
reasonable if the interaction between the subsystems is small, i.e. if their lengths are 
large compared to the length scale associated with localization in the linear problem: 
Lq £ (£ is the localization length). On these scales the various subsystems are 
statistically independent. This is the content of (|35|) . It motivates the definition of 
the L- independent quantity, 

i?(p,W0 = [P( P ,W r ,L)] 1 / i . (36) 

This scaling relation was verified numerically, starting with a uniform distribution 
and using periodic boundary conditions. It was found that for lattices of sizes 
16 < L < 128, R is independent of L. Therefore, to calculate P it is sufficient to 
evaluate, 

P (p,W) = P(p,W 1 L ), (37) 

for a system of some size Lo- It is assumed that the scaling found, will hold 
for systems of arbitrarily large size, which allows to increase also the localization 
length. It is convenient to perform a transformation to a new quantity Q(p, W) as 
Q = P / (1 - P ), which yields, 

1 
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Figure 3. (Color online) The function Q(p, W) (top) and the same function in 
scaled coordinates (bottom). The black dashed and red dotted lines, are showing 
asymptotics for small and large arguments of q, have slopes £ = 9/4 and r) = 5.2. 
(Fig. 4 and 5 of [54]) 



The limits p — > (linear problem) and W — > (non-random system) are singular, 
therefore it is expected and indeed observed that the function Q(p, W) is not an 
arbitrary function of p and W, but it can be written in a scaling form 



w aq (~> 



where q(x) is a singular function at its limits, q(x) ~ cix~^ for small x. while 
q(x) ~ C2X~ V for large x. The top of Fig. [3] collapses to one curve as shown in 
the bottom of Fig. [3] This is the numerical justification for (|39|) . It also provides the 
values of the exponents a = ot\ « 1.75, C ~ 9/4 = 2.25, r\ w 5.2, c\ « 2.5 • 10~ 7 and 
c 2 w 1.8- 10- 18 . 

In particular, for a fixed density, P c h ~ /0 9 ^ 4 , since the probability of chaos is 
fixed and non-zero, spreading is expected. In this aspect model (TTJ) differs from the 
corresponding N— body problem, where for a fixed density and in the thermodynamic 
limit, localization was found |35l I36| . Now let us assume that we consider the states 
with the same fixed norm Af on lattices of different length L. Then, p = Af/L, and 
one finds 

Xi-CyyCMMi-C) _ L-5/W9/4 

This quantity, as expected, grows with the norm Af and decreases with the disorder W. 
We see that because £ > 1, the probability to observe chaos in large lattices at fixed 
norm tends to zero. This result may have implications for the problem of spreading of 
an initially localized wave packet in large lattices. In this setup the norm of the field 
is conserved, and the effective density decreases in the course of the spreading. This 
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means that as a function of time L increases and P c h decreases, therefore, spreading 
takes place for less and less realizations of the disorder. 



2.3. Study of the relation to phase-space structures 

It is expected that in the course of spreading, as the nonlinearity decreases, the 
trajectories in the ip (x), if)* (x) phase-space will become more and more regular. In 
particular, trajectories that look chaotic on some scale may eventually look regular. 
For this purpose the technique of the time-dependent Lyapunov exponents was 
introduced [55] . It suggests that an initially chaotic wavepacket may stick to KAM like 
trajectories resulting in a slow-down of the spreading. A mechanism for spreading that 
involved a resonance of three oscillators resulting in a mechanism for Arnold diffusion 
was recently proposed [56 1. It is a spot in phase-space where the local Lyapunov 
exponents are positive and there is wandering in phase-space, presumably, in the 
regions where tori are destroyed. 

A recent qualitative argument for a finite probability for spreading, combined 
with numerical results was presented [57 . It is instructive to compare it with rigorous 
statements on the same model [39] . 



2.4- The renormalized perturbation theory 

Since there is no spectral theory for nonlinear equations analysis in the framework 
of the time dependent perturbation theory was performed [53 031 Q3] . The objective 
is to develop a perturbation expansion of the c m (t) of in powers of (3 and to 
calculate them order by order in j3. The required expansion is 

c n (t) = 4°) + /3c« + /3 2 c( 2) + ' • ' + P N ci N] + Q n , (41) 

where the expansion is till order N and Q n is the remainder term (clearly, ch, and Q n 
are random variables) . The initial condition 

c n (t = 0) = 8 n0 , (42) 

was assumed. The equations for the two leading orders are presented in what follows. 
The leading order is 

= S n0 . (43) 
And the first order is 

/ 1 p i(E n -E„)t\ 

^-^(hbir-)- ,44) 

The divergence of this expansion for any value of /? may result from three major 
problems: the secular terms problem, the entropy problem (i.e., factorial proliferation 
of terms), and the small denominators problem. To eliminate the secular terms the 
ansatz (TTU1) is replaced by [TJ] 

i/> (as, t) = Y,cn (t) e- iE '^u n (x) , (45) 

n 

where 

E' n ^E^+f5E^+^E^+.-. (46) 

and E { n 0) are the eigenvalues of Hq (clearly, En\ are random variables), and E' n are 
the renormalized energies. The new equation for the c„ is given by 
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id tCn = - E' n ) c n (47) 

Inserting expansions (|4"Tj) and (|4"rH) into (14"T1) and comparing the powers of (3 without 
expanding the exponent in (3, produces the following equation for the k — th order 



E ^« mir 



mi 7712^3 



fe— 1 fe— 1— r fe— 1 — r— s 

/ > / y / j Sii 

r=0 s=0 i=0 



The relation to the order by order expansion in powers of /? was presented in [T3] . 

The perturbative approach gives an explicit value of the wavepacket for long times, 
with a rigorous control of the error (Fig. [3]). In particular, it rules out spreading for 
the NLSE with J = 1, w = 4, /3 = 0.1 for time of order t = 10 5 [5S]. 

2.4-1- The remainder of the expansion In order to control the solution one has to 
control, Q n , the remainder of the expansion (|41l) that can be written in the form 

Cn (t) = C„ + Q n , (49) 

with 

Cn=Y.^ C n- ( 5 °) 

The linear part of the full equation for the remainder is given by 

id t Ql U = W n (t) + E M nm (t) Q 1 ^ + E Mnrn (t) * . (51) 

m m 

The contribution of the remainder term is 

|Q n | < const • e ^+NlnP+lnt e -^ n \^ ^ 

where 7 = 1/^ is the inverse localization length, x n is the localization center 
of the n— th state. Note that for a given t and /3 there is an optimum N for 
which the remainder is minimal. Additionally, for any fixed time and order N, 
lim^o \Qn\ I P N ~ 1 = (see definition at [13]), which shows that the series is in fact an 
asymptotic one |60| . Using a bootstrap argument, one can show |13) that until some 
time t* , the dynamics is governed by the linear part and the remainder is bounded 
by 

\Q n (t)\ <A(3 N t-e-^, (53) 

where 7 is the inverse localization length and A is a constant. Therefore to estimate 
the remainder one can integrate (1511) at least up to t*. It is useful to integrate up to 
some large time, t t* , and then to extrapolate using the linear bound (|53[) up to t* . 
In the next subsection it will be proposed how to determine t* in practice. 

For time shorter than t* there is a front x (t) oc In t such that for x n > x (t) both 
the remainder, Q n (t) and c n (t) are exponentially small. 
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Figure 4. log 10 t„ as a function of /3 1 for different orders. 4th order (solid 
blue), 3rd order (dashed green) and 2nd order (dotted red). The parameters are: 
w = 4, J= 1. 



2-4-2. Practical estimate of the remainder In this section it will be demonstrated, 
how the numerical scheme for calculations in the framework of the perturbation theory, 
is implemented in practice for a specific realization of the random potentials {e x }. It 
was found that there are some modes that result in the largest contributions to the 
remainder. These modes are located near the initial mode and therefore can be easily 
identified and subtracted. For numerical evaluations one defines a norm 



the definition of ||Qzm|| 2 is similar. It was found that for ||Qzi ra || 2 < 0.1, ||Qiin|la ^ s 
close to ||Q || 2 ) which leads to the definition of t*, 



For small nonlinearity strength, ft, i* is very large and therefore the integration of ([BTj) 
to t* is very time consuming. Equation (|55[) can be used to extrapolate linearly from 
the time interval where (|5 1 [) is solved to i*. Practically, one can find t* from ([55)1 by 
extrapolation. In Fig. 3] a plot of log 10 i* as a function of ft^ 1 for different orders is 
presented. 

A systematic improvement with the order of the perturbation theory is found. 
Note, that even with moderate nonlinearity strengths, namely, ft < 0.08, one can 
achieve a good approximation of the solution up to very large times. In this way the 
perturbation theory combined with the solution of the linear equation (|5ip and the 
criterion (|55p may be used to obtain the solution of the original equation ([1]) up to 
t < t*. For small ft the time i* is very long, as is clear from Fig. 0] Removing the 
dominant modes, which were mentioned in the beginning of this subsection, allows to 
obtain reliable results for times larger by more than one order of magnitude from the 
times presented at Fig. 0](see, Fig. 10 in [H]). When one considers the smallness of ft 
one should consider actually the smallness of fte 2 due to the exponential proliferation 
of the number of terms (see Eq. (4.6) of |13|). 




(54) 



\Qun (*,)|| 2 =0.1. 



(55) 
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3. Rigorous results 

The rigorous analysis of the NLSE equation with a random potential turned out to 
be very difficult as well. The results so far are very limited in scope, yet of sufficient 
interest to point out the nature of the problem at hand. Most notably we have the 
following conclusive results. 

3.1. Many body localization 

Consider the system of N interacting particles via finite range interactions, on a finite 
d— dimensional lattice, Z d , with the Hamiltonian 

N 

(w) = ]T [-A i + XV ( Xj ;u)] + J2 Va ( Xi - x ) (56) 

3=1 i<3 

acting on the Hilbert space Ti^ N ' ) = £ 2 ('E d ) N . It is assumed, for simplicity, that the 
Vij are all uniformly bounded by a constant, < a < oo, and the random potential 
V (xj)ui) is given in terms of a collection of i.i.d random variables, {V (x,w)} ieZ d, 
with bounded probability distribution and finite moment generating function 

E(exp(tV (0))) < oo V*. (57) 

Under the above assumptions on the potentials, it was proved that for a large 
disorder, the spectrum of JjW is a dense point spectrum with probability 1, and 
the eigenfunctions are exponentially localized, after an appropriate distance function 
is defined on the lattice, between clusters of particles 138,125]. This important result 
shows that the inter-particle interactions between, say, electrons in a solid, cannot 
destroy Anderson localization, at least under some favorable situations. However, 
this result requires the size of the disorder to be dependent on N, the number of 
particles! Therefore, it cannot be applied to systems with non-zero density of particles, 
as systems studied in |35[ |3T>] . 

3.2. Quasi-periodic perturbations 

We now turn to the fully nonlinear problem. The first important result in this 
case is |47| . where the time- independent problem was explored. The nonlinear term 
of the NLSE is considered as a small perturbation of the linear dynamical system 
corresponding to the linear Anderson problem on the lattice. One is then led to 
consider the KAM theory in infinite dimensional phase-space, as a way to construct 
periodic and quasi-periodic solutions (in time) for such models. This has been shown 
to apply to models with good Diophantine properties (linear combination with integer 
coefficients, bounded away from zero) of the eigenvalues of the linear part. The 
possibility that the solutions are localized in space and quasi-periodic in time, notably 
leads to the question of whether linearizing the NLSE around such a solution, will 
result in a linear system with localized states only. The idea behind the use of quasi- 
periodic in time models is that it comes from a formal iterative scheme for solving the 
NLSE with localized solutions. If a localized solution of the equation is assumed, it 
can be expanded in terms of the normal modes (eigenfunctions) of the linear problem, 
so that 

b(N) 

^ N ( x ) = ^ CjUj (x) e~ iE i\ (58) 

3=1 
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where N is the iteration number, and b (N) is the number of modes generated by N 
iterations. In this approximation the nonlinearity is given by, 

2 



A |^jv| 2 = A 



b(N) 
3 = 1 



(59) 



Therefore, the approximate dynamics is governed, to this order by, 

id t ip= (-A + VOV + AI^M 2 ^. ( 60 ) 
A uniform (in TV) proof of complete dynamical localization, for the linear equation 
presented above, would imply localization for the NLSE problem. However, in this 
iteration scheme b (N) grows very fast in N, and furthermore, the frequencies in 
A|?/>Ar| may be of an arbitrarily small value, which is hard to control |61| . This is 
just another manifestation of the small denominator problem one encounters in other 
perturbative methods. 

This approach was followed in a series of papers, beginning with [52], where the 
Anderson model on a lattice is perturbed by an exponentially localized time-periodic 
potential. It was shown that the corresponding Floquet operator has purely dense 
point spectrum, with exponentially localized eigenfunctions. This implies, by general 
spectral theory that the original equation has dynamical localization, namely, the 
initial solution does not spread. This was further generalized to the much harder 
case of exponentially localized, quasi-periodic in time potential perturbations of the 
Anderson model on the lattice with a similar result - complete localization for small 
potential perturbations |61| . The key estimates can be formulated as localization for 
the Floquet operator of the type (for two non-commensurate frequencies) 

d d 

K = -ioj x — — iw2T^-+eA+V r +T^i(j)cos27r0i-(-VT2(j)cos27re'2,(61) 

Oui UU2 

acting on the Hilbert space Hk = (g> L 2 (T 2 ) (here, T 2 stands for a two- 

dimensional torus). Then localization was proved with large probability, and for a 
corresponding set of Diophantine conditions needed in the quasi-periodic case [SU [HI] . 
In particular, it was shown that with large probability, K has a dense point spectrum, 
with exponentially decaying eigenfunctions. However, the estimates deteriorate as 
the quasi-periodic frequencies approach zero, which makes them hard to use for the 
solution of the full NLSE problem. 

Further results on models approximating in various ways the NLSE problem, 
include [12] > where it is shown that quasi-periodic solutions in time exist, and that 
they bifurcate from the corresponding quasi-periodic solutions of the linear equation. 
Another result is that if the nonlinearity coupling constant vanishes at space infinity 
(at some polynomial rate), then the NLSE with Anderson potential on the lattice has 
only localized solutions for large disorder [7]. 



3.3. NLSE with a random potential 

Turning our attention to the full NLSE with a random potential problem, very little 
is known. The most important result, which indicates localization for large times with 
large disorder and small nonlinearity is in the work of [9]. The main result can be 
formulated as follows. The small parameter in the problem is e = J + j3, and w = 1. 
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E 



Therefore, small e implies strong disorder (J <C 1) and weak nonlinearity (/3 <C 1). For 
an initial data in £ 2 , localized at the origin, in the sense that, 

\^(j,t = 0)\ 2 <8<1, (62) 

\]\>R 

it was proven[3J, that for all A and 6 > there exists a constant C (A) > 0, e (A) > 
and K {A) > A 2 such that for all t < (S/C) e~ A , 

E hW)| 2 <2£, (63) 

\j\>R+K 

with probability 

l-exp(-|e- 2 < /CA ), (64) 

for all < e < e (A). In this context, it should be noted that the methods used to 
prove such results, follow the infinite dimensional generalization of dynamical systems 
theory. It includes a construction of normal form transformations to approximately 
"diagonalize" the system. Such methods were used already in [49] to prove a similar 
result for the nonlinearly coupled random masses. This result implies that the 
£ 2 — norm grows at most logarithmically in time [pj. However, it is hard to compare 
these estimates to other results, since the constants are not controlled in this analysis. 
Note, that the corresponding results (|52|) and (|53|) of the perturbation theory are for 
arbitrary strength of disorder, and improve when it increases. 

Another class of results, which give explicit estimates on the solution of NLSE with 
a random potential was developed in [58, I13[[T4"] . One can construct renormalized time- 
dependent perturbation theory for the solutions, starting from the eigenfunctions of 
the linear problem as a basis. It is then shown, that after eliminating the secular terms, 
we get an explicit expansion, which is computable. Each term depends explicitly on the 
potential of the linear problem and the eigenvalues. Therefore, by using the known 
and some new results on the linear spectral problem it is possible to estimate the 
various terms in the expansion on average. In particular, in |63| it is shown that for 
the linear problem in one dimension, if the potentials are uniformly bounded by some 
finite constant, then for all potentials the minimal distance between the eigenvalues 
is nonzero, bounded below by e~ cN , where N is the size of the lattice. This is then 
used to control rigorously the first order term in the above expansion (|4]~|) . (|44|) . and to 
conclude that to this order the exponential localization persists in the nonlinear case. 
It should be noted that the lower bound on the distance between the eigenvalues, 
though exponentially small, is sufficient to the purpose of controlling the various 
terms in the renormalized expansion. It also allows to obtain bounds without the 
usual Diophantine conditions. Moreover, it implies a bound on the expectation of the 
derivative of the eigenfunctions with respect to the potential, 

dtl)i (x) 



E 



E 
l 



de 3 



< CN. (65) 



However, to control the higher order terms in the expansion, a control of linear 
combinations of more than two eigenvalues is needed. In the work of |34J there are 
Diophantine estimates on the eigenvalues of the linear problem that can similarly 
control many other terms in the expansion of [581 [TBI IT4"] . However, this is not 
sufficient for complete control in the probabilistic sense of the linear combinations 
of the eigenvalues of the linear problem, to bound all the terms in the expansion. 
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4. Numerical results 

The main problem one encounters with the numerical calculation for chaotic systems is 
the exponential sensitivity to numerical errors. For this reason one cannot perform the 
calculation for substantially long time scale with a control on the errors of a specific 
solution. The validity of the results is traditionally tested by changing the size of the 
steps verifying that the results are not affected. The correct result is assumed to be 
the one that is found in the limit of a vanishing step size. The main problem with such 
assumptions is that the limit of zero time step may be singular. Moreover, there is 
no theoretical understanding that the numerical solutions are some types of sampling 
of the phase-space as is the situation for chaotic systems. Typically, one finds that 
after some time, spreading starts, as shown in Fig. [2] Eventually, it turns into a 
sub-diffusion and the second moment grows according to (|26| with the exponent l /s 
[101 [HI [T2 , 64. 65j. The wavepacket is typically of the shape presented in Fig. [TJ A 
theory similar to the one presented in Subsection 12.11 but with the nonlinear term (j4]) 
was developed and was found to agree with the numerical results [SJ Q21 UHl HH1 ES] • 
The longest time for which numerical calculations were performed is t — 10 8 (in units 
where J = 1 in ([T])). 

In order to follow the dynamics of a wavepacket, typically the SABA algorithm is 
used. This algorithm belongs to the family of split step algorithms and evaluates the 
wave packet in small steps, changing from coordinate space to momentum space. The 
disorder and nonlinear interaction are applied in the coordinate space, the wave is then 
transformed to the momentum space, where the kinetic energy term is applied. The 
solution in real space is recovered by transforming it back to the coordinate space, and 
the procedure is repeated. Nearly all numerical calculations for this problem use such 
methods. Additional details on the SABA algorithm, can be found in [12J. Like any 
numerical algorithm, the SABA algorithm accumulates errors during the calculation 
which grow with the time of the integration. There is no reliable estimate of these 
errors. 

Double humped states were studied numerically in the presence of nonlinearity 
that is not too strong. It was found that the spreading of a wave packet prepared 
initially near some site O is substantially stronger if there is a double humped state 
with one of its humps near O, than if the states peaked near O are single humped. 
It was found |67) that there is a regime of parameters where /3 is sufficiently small so 
that the double humped structure is preserved but the packet is not only oscillating 
between the humps but also leaks to other states, leading to spreading. Additionally, if 
j3 is small enough so that the oscillations between the two states are not suppressed in 
the double- well model, then the double humped states will contribute to the spreading 
for the NLSE. But since double humped states are suppressed and do not contribute 
to the spreading for high nonlinearities, it cannot be claimed that they dominate the 
spreading for the NLSE. 

5. Discussion and open problems 

All numerical calculations exhibit spreading of an initially localized wavepacket 
(although some part of it may not spread). The spreading results in sub-diffusion 
with the exponent All rigorous and analytical theories predict that asymptotically 
the spreading cannot be faster than logarithmic in time. The main difficulty is, that 
there is no regime of parameters, where analytical and numerical results agree for a 
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long time. For short times (t < 10 4 ), perturbation theory was found to agree with the 
numerical results (and for t < 10 5 in |59|). 

5.1. List of open problems 

We list the questions that may be explored by the various communities: 

(i) What is the asymptotic time scale? Namely, what is the time scale where the 
theories predicting suppression of sub-diffusion become effective. 

(ii) When does spreading start and how this time depends on parameters? 

(iii) Is it possible to analytically derive the scaling theory presented in Subsection l2.2l .'' 

(iv) Can one prove the bound of the terms of the perturbation theory presented in 
Subsection 12.41 ? In this theory it is difficult to control denominators of the form 
/; = Y]j—i CiEi where Ei are the eigenvalues of the Anderson Hamiltonian, Hq, 
and Ci are integer constants. Numerical calculations show that these satisfy a 
central limit theorem. In particular, one would like to prove that, 

(|7t) <0 ° for0 < s < 1 - (66) 

(v) Is the system chaotic? Namely, is there an exponential instability of the motion 
in the ip (x), ip* (x) phase-space? This is a fundamental question since this phase- 
space is of infinite dimension. 

(vi) Are there KAM tori? Is there sticking to these tori and is it stable to the numerical 
errors? 

(vii) Can one design experiments that can be extended to the asymptotic long time 
regime? 

(viii) May the control of the numerical scheme be improved? 

(ix) Another possible mechanism for spreading is tunneling to exponentially large 
distances with exponentially small probability in space and time. This cannot 
be decided using numerical calculations, since only relatively small systems are 
accessible. Another method should be developed to resolve this issue. 
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